Hitherto, in system identification based on a parametric model, an objective function for parameter estimation including a prediction error is constituted, and a process of deciding parameters to minimize the objective function is performed. The objective function that is most normally used is a quadratic function. In this case, the process of deciding parameters to minimize the objective function is called “least-squares method estimation.”
Parameters to be identified are defined as x=(x1, . . . , xm) and plant models are defined as {circumflex over (ƒ)}(x).
A problem of deciding parameters using the least-squares method is defined by the following expression as a problem of deciding parameters which allow square error ε2 with respect to observation values ƒi (i=1, . . . , n) to be minimized. Here, m and n are natural numbers and m<n is satisfied:
      ɛ    2    =            ∑              i        =        1            n        ⁢                  (                              f            i                    -                                    f              ^                        ⁡                          (              x              )                                      )            2      
Then, it is a problem of identifying m parameters on the basis of the following m simultaneous equations obtained by performing partial differentiation on the square error ε2 by using the individual parameters. The simultaneous equations are defined as normal equations where
            ∂              ɛ        2                    ∂      x        =  0.
Thus, the problem to be solved by using the least-squares method is defined as follows:
                    [                              ∑            n                                                          ⁢                                    (                              -                                                      ∂                    ɛ                                                        ∂                    x                                                              )                        ⁢                          (                              -                                                      ∂                    ɛ                                                        ∂                                          x                      T                                                                                  )                                      ]            ·      Δ        ⁢                  ⁢    x    =            ∑      n                            ⁢                  (                  -                                    ∂              ɛ                                      ∂              x                                      )            ⁢              ɛ        .            
Here, the following expression is obtained when terms of partial differentiation are expressed by matrix A:ATA·Δx=ATε.
When this expression is transformed by replacing Δx by {circumflex over (θ)}, the following expression is obtained:{circumflex over (θ)}=(ATA)−1ATεwhere ε is the difference between observation value ƒi and plant model {circumflex over (ƒ)}(x) and {circumflex over (θ)} is the estimated value of the parameter to be identified.
Particularly, in the least-squares method, it is known that the term ATA of the following matrix is a sparse matrix. Here, it is assumed that matrix A is an m×m matrix, aij represents an ij element of matrix A, and n represents the number of measurements.
            A      T        ⁢    A    =      [                                                      ∑                              b                =                1                            n                        ⁢                                          a                                  b                  ⁢                                                                          ⁢                  1                                            ⁢                              a                                  b                  ⁢                                                                          ⁢                  1                                                                              …                                                    ∑                              b                =                1                            n                        ⁢                                          a                                  b                  ⁢                                                                          ⁢                  1                                            ⁢                              a                bm                                                                          ⋮                          ⋱                          ⋮                                                                ∑                              b                =                1                            n                        ⁢                                          a                bm                            ⁢                              a                                  b                  ⁢                                                                          ⁢                  1                                                                              …                                                    ∑                              b                =                1                            n                        ⁢                                          a                bm                            ⁢                              a                bm                                                          ]  
Generally, a person who performs modeling is different from a person who writes analysis code. Thus, studies of how to efficiently perform calculation when sparse matrix ATA is given have been made. However, the form of a sparse matrix significantly varies or properties thereof locally vary in individual object problems, and thus it is necessary to perform a process of determining a storage scheme by testing some storage schemes at the time of calculation or to perform an experiment for determining a storage scheme in advance. Also, since the same storage scheme is selected for the entire sparse matrix ATA, it is difficult to efficiently store the sparse matrix.